New approach to calculating the News 
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We consider the problem of finding the gravitational radiation output, or news, within the context 
of a numerical simulation of a spacetime by means of the null-cone, or characteristic, approach to 
numerical relativity. We develop a method for computing the news that uses an explicit coordinate 
transformation to a coordinate system that satisfies the Bondi conditions. The method has been 
implemented computationally. We present results of applying the method to certain test problems, 
demonstrating second order convergence of the news to the analytic value. 
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I. INTRODUCTION 

In many simulations in numerical relativity, an impor- 
tant objective is to compute the emitted gravitational 
radiation. Within ADM- type approaches, one uses com- 
puted data near the outer boundary of the simulation to 
estimate the radiation at future null infinity in the 
form of an expansion of spherical harmonics. The un- 
derlying theory was developed some time ago Q, 0, > 
and the computational implementation is discussed, for 
example, in |j. 0.Q However, this paper uses the char- 
acteristic, or null-cone, approach with Bondi-Sachs coor- 
dinates, in which T + is included in the compactified grid, 
and therefore the gravitational radiation, or news, can in 
principle be found exactly |2j. 

Even so, the computation of the news is not a straight- 
forward matter. The news takes a very simple form, 
N = ij r £ u , when the null-cone coordinate system sat- 
isfies the Bondi conditions (the meaning of these terms 
and symbols will be given later). However, the Bondi 
conditions are conditions at X + (where I = 0) , and in any 
physical simulation the freedoms in the null-cone coordi- 
nates are set in the interior of the spacetime. In general, 
a computed metric will not satisfy the Bondi conditions. 
Thus, a procedure for computing the news in a general 
Bondi-Sachs coordinate system has been developed Q. 

A motivation for this work is that some difficulties 
with the news computation of have been reported , 
specifically significant gravitational radiation when the 
gravitational field is static. The cause of the difficulties 
is not known, and it may be that they can be resolved by 
technical improvements to the method and code |lfj| , but 
as yet no such results have been reported. Our method 
has also not been tested for the problem of 9] , and so it 
is not certain that it will be of assistance; but because 
of its simplicity (see below) it would be much easier to 
analyze, and thereby resolve, any problems. 

However, our primary motivation is that the method 
presented here is much simpler than the previous method. 
Conceptually, the idea is straightforward: we construct 
a transformation to the "right" coordinate system, and 
consequently there is a means to monitor the accuracy 
of results obtained. The technical differences between 



the two methods are summarized in Sec. IVII Here we 
note that the present method requires two fewer evolu- 
tion equations than the previous method, and also that 
the most complicated formula in the present method is 
Eq. (|25|l . whereas the formulas for the previous method 
were given in an appendix taking up a whole page. In 
numerical work, simplicity means that there is less op- 
portunity for the introduction of error. Of course, the 
previous method has been tested against available ana- 
lytic solutions (as will be done for the method presented 
here), but these solutions are rather special. Thus, what- 
ever the computational performance of the two methods, 
it would still be necessary to develop the present method 
and use it to validate results obtained by the previous 
method. 

We construct an explicit coordinate transformation be- 
tween the general Bondi-Sachs coordinate system and a 
Bondi-Sachs coordinate system that satisfies the Bondi 
conditions (given in Eq. (JjJJ). The coordinate transfor- 
mation on X + is known |8( , but constructing the transfor- 
mation throughout the spacetime is much more difficult. 
The key point of our method is that, while it is in prin- 
ciple possible to find the transformation throughout the 
spacetime, we do not need it. Rather, we need to eval- 
uate \J ,iu on I + (where £ = 0), and thus we need the 
transformed metric off T + only to first order in t. We 
proceed by expressing the coordinate transformation as 
a Taylor series in £. The algebraic calculations are quite 
lengthy and require the use of computer algebra. 

Our news calculation has been implemented computa- 
tionally, and used in problems in which an analytic result 
is known. Results are presented below, and it is seen that 
the errors are second-order convergent to zero. 

The paper starts with a summary of relevant results 
and notation for the characteristic formulation of numer- 
ical relativity (Sec.[Hj|. The coordinate transformation 
and computer algebra results are given in Sec. IIIII and 
then the procedure for computing the news - at both an- 
alytic and computational levels - is described in Sec. II VI 
The computational tests and results are presented in 
Sec. The Conclusion, Sec. IVII gives a comparison 
between the method in this paper and that of and 
also discusses possibilities for further work. 



2 



II. NOTATION 

The formalism for the numerical evolution of Einstein's 
equations , in null cone coordinates, is well known 0, |1| 
(see also [Tl|, [T^, HH, |l4|). For the sake of completeness, 
we give a summary of those aspects of the formalism 
that will be used here. We start with coordinates based 
upon a family of outgoing null hypersurfaces. We let u 
label these hypersurfaces, x A (A — 2,3), label the null 
rays and r be a surface area coordinate. In the resulting 
x a = (u, r, x A ) coordinates, the metric takes the Bondi- 
Sachs form [HE] 



ds 2 



W 



e^(l + — ) - r 2 h AB U A U B du 



-2e 2p dudr - 2r A h AB U a dudx A + r 2 h AB dx A dx a , (1) 

where h AB h B c = S A and det(h AB ) = det(q A B), with 
q A B a unit sphere metric. We work in stereographic co- 
ordinates x A = (q,p) for which the unit sphere metric is 



q AB dx A dx B 



-^{dq 2 +dp% 



where 

P=l + q 2 +p 2 . 
We also introduce a complex dyad q A defined by 



P 



(l,z), q A = —(!,«) 



(2) 



(3) 



(4) 



with i = \/— T. For an arbitrary Bondi-Sachs metric, h AB 
can then be represented by its dyad component 



J = h AB q A q B /2, 



(5) 



with the spherically symmetric case characterized by 
J = 0. The full nonlinear h AB is uniquely determined 
by J, since the determinant condition implies that the 
remaining dyad component 



satisfies 1 
field 



K 2 



K = h AB q A q B /2 (6) 
J J. We introduce the spin-weighted 



U = U q A , 



(7) 



as well as the (complex differential) eth operators 8 and 
o see 16] for full details). 

The news calculation is performed in conformally com- 
pactified coordinates. Specifically, (u,r, x A ) — + (u,£, x A ) 
where I = 1/r. In (u,£, x A ) coordinates, the compacti- 
fied metric is ds 2 = £ 2 ds 2 where I is a conformal factor 
with future null infinity X + given by £ = 0. The com- 
pactified metric will be denoted by g a @ , and the general 
compactified Bondi-Sachs metric is 



,-2/9 



^10 



-2/3 



9 AB 



Va, 9 
AB 



1A 



-OA 



-2PttA 

a 



r 

^00 _ 



0, 



(8) 



where V a = £ 2 (l + £W). In addition to the general com- 
pactified Bondi-Sachs metric Eq. ©, we will refer to 
the compactified Bondi-Sachs metric satisfying the Bondi 
conditions, and such quantities will be denoted with a 
suffix , with compactified metric and coordinates g?™ 
and (u[ B ],£[ B ],x A B ^), respectively. OnI + , i.e. £[ B ] = 0, 
satisfies the Bondi conditions 



9[ B] = o, gfi } = 0, 9° [B] = 1, 9$ 



n AB 
1[B] > 



(9) 



where q A ^ is a unit sphere metric with respect to the 
Bondi angular coordinates x A B y 



III. COORDINATE TRANSFORMATION 



We define a coordinate transformation, near I + , be- 

A 1 



tween (u 7 £,x A ) and £[ B ], x \b])'- 



i[B] 



i [B]Q + A u £ + C u £ 2 , £ [B] =lu£ + C l £ 2 



where 



w, x 



B10' 



A* 



l>[B]Q> 



xf B]0 + A A £ + C A £ 2 , (10) 



A u are all functions of 



and u only. It will turn out that the C a are irrelevant 
to the news calculation, but they are included so that 
the Jacobian (constructed by differentiating the above) 
is manifestly correct to first order in I. We also introduce 
complex quantities 



A = q A A A and X 



Qax a B ] Qi 



(11) 



both of which are defined to have spin weight 1. The 
general metric quantities are expanded as a Taylor series 
to first order in £, with the o quantities being the values 
at £ = 0, i.e. on X + , and the e quantities being the first 
I derivative evaluated at £ = 0, i.e. (di)e—o. 

= 0o U = U + £U e , J = J + £J e , 

K = K +£Ki, V a =W at . (12) 

The same first-order expansion in £ and notation will 
be used for the metric quantities in the coordinate sys- 
tem satisfying the Bondi conditions (e.g., Jt B i = J[ B ]o + 
£J[B]e)- The general compactified Bondi-Sachs metric, 
and the Bondi-Sachs metric satisfying the Bondi condi- 
tions, are related by 



~af3 
■9[B10 



dx" B ] dx 







(13) 



with the factor uj 2 appearing because there is also an 
implicit change of compactification factor from £ to £\m . 

We have used Maple to find to first order in £, and 
in doing so we have imposed the conditions (c.f. Eqs. 
(39) and (40) in @). 



(d u + U B d B )x A ]0 = 



(14) 



3 



and 

(d u + U^d B )u [B]0 ^ioe 2(3 °. 
We found that the Bondi conditions 



9[B]o — 0> 9[b]q — 



i A 



?[B]0 



= 1 



(15) 



(16) 



are satisfied identically. The remaining Bondi condition 
is 



p.AB _ AB 
9[B]0 — 1[B] ■ 



(17) 



Taking the determinant of Eq. (|17Jl leads to an explicit 
expression for uj 



Pu 



<l[B},qP[B}, P - q[B],pP[I 



(18) 



where P[ B ] = 1 + <7[b] 2 +P[b] 2 - The remaining content of 
Eq. 117|) can be written as 

J[B]o = -ff[sf 9A[B]9s[s]/2 = 0, (19) 

where qa[b] is the complex dyad appropriate to the met- 
ric of^- The Maple calculation finds 



BIO 



AP 2 



U X 2 C 2 - 2K dXdX + 4J XBX( 

-4K XQX(+ J (dX) 2 + J (BX) 2 
+ P 2 XX) 2 lu 2 ] = 0, 



(20) 



where Q = q — ip. Eq. I|2U[) is not used to determine any of 
the transformation parameters, but rather as a check on 
the accuracy of the computation, as described in Sec. lIVI 



In addition to satisfying the Bondi conditions, g?^ is 
also a Bondi-Sachs metric and must satisfy the conditions 



9[B] 



= 0, 



-0.4 
9[B] 



= 0. 



(21) 



Now, <?[^] = and g?^ = can be expressed as linear 
equations in A u and C u respectively, but these values 
will not be needed. Also, gP™ = and g^ e — can be 

expressed as linear equations in A A and C A respectively: 
the expression for C A will not be needed, but the expres- 
sion for A A will be used in the news calculation and, in 
complex form, 

A = [J 8u [B] (2X( + dX) + J du [B] BX 
-K Q {nu [B] (2XC + BX) + Bu [B] dX)}/(2uj). (22) 

The computer algebra calculation evaluates the inverse 
Jacobian at £ = 0, from which we find a result that will 
be used in the news calculation 



(u,£,x A ) = ^—(l,0,U A ). (23) 



du 



IB] 



The computer algebra yields an expression for J[ B ^f. 
Defining 



Tx = J 3A - K dX, 
T 2 = J Q {BX) 2 + dX(J dX ~ 2K BX), 



(24) 



then 



J [B]e = P 2 (4(l - P [B ])AT 2 + 2XP [B] (J e (dX) 2 + J e {dX) 2 + 2(MTx + 2(XJ M + dA(J 3X - K BX) - 2(XK dA 
-K$X(2(X + BX)) + 4(J e X((X + QX)) - AAX{-4( 2 J X + P [B] (e -2 * (2(U e X + (U e BX + U$X) + 2A. U + 
AV a e + (U BA + U SA) + 2(AU ) + 2fT x )) - P 2 X 2 A(4(X((J X + Tx) + T 2 ) + l&CAXTx)/{%XP [B fu 2 ). (25) 

I 



IV. PROCEDURE FOR COMPUTING THE 
NEWS 



We can now describe the procedure used for calculating 
the news. First, we look at the process from an analytic 
viewpoint, and then we give some details of the computa- 
tional implementation. As discussed in Sec^ we assume 
that the general Bondi-Sachs metric is known throughout 
the spacetime. 



A. Analytic procedure 

1. Solve Eq. lO to find x A B]Q = x A B]Q (u,x A ). 

2. Eq. (HHJ) gives u = lj(u,x a ). 

3. Use Eq. I|2UI) to find by how much J\ B ]o differs from 
zero, and thus obtain an estimate of the accuracy 
of the calculation. 

4. Solve Eq. H15|) to find M[b]o = U[B]o( u i x )■ 

5. Eq. (EH) gives A = A(u, x A ). 
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6. We note that 
/ dJ\ B 



\d*[B] 



dx 11 8J[b] 



! [B ]=0 

d£ dJ<B 



d£ [B ] dx» 

J[B]t 



di [B] di 



to 



(26) 



because Jrgi is zero on T + , and so derivatives in 
directions other than I vanish. Then the news N = 
N(u,x A ) is 



N 



1 d 2 J [B] 



[B] 



2 dl[B\du 
f du d 
\du[ B ] du du[ B ] 



dx 1 



d N 


\ J [B]t 


du [B ] , 


1 2uj 


■ 8 ) 


J[B]t 


dx A J 


2uo 



(27) 



The terms du/duis] etc. are evaluated at i = 
and were given in Eq. above. 

7. Finally, we need to express the news as a function 
of Bondi coordinates, i.e. to transform N(u, x A ) to 
N(u [B] ,xf B X 



B. Computational implementation 

The news module has been written to interface directly 
with the null gravity code in its current form |17|. Thus 
we use a compactified radial coordinate x = r / (R + r) , 
with R = 1. There are n x points in the x direction in 
the range [0.5, 1] (corresponding to 1 < r < oo). We 
use stereographic angular coordinates, with the number 
of grid points in each angular direction n n . 

Instead of solving equations l|14l) and l|15|) . we convert 
them to ordinary differential equations as in reference § , 
and solve them by means of the second order method 



(x A(n+1 *> - x A ^)i = Au(U A ot +k2) , 



(n+l) 



L (n) ) 
[B] ' 



Au (toe 



(28) 
(29) 



At the initial time-slice we match the Bondi coordinates 
to the Bondi-Sachs coordinates: i.e. at u = 0, U[ B ] = 

J i 



and xA, 



This also gives A = and Ji 



[B]t 



on the initial surface. The solution of Eq. (|28|l leads to 
x A = x A (u, x A B ^), and this needs to be inverted to give 



[b] 



:j]jj(u,a; ). At this stage the inversion has not 
been implemented in general, because it is not needed in 
the test examples considered below. 

Second derivatives do not appear directly in the formu- 
las leading to the news. First angular (eth) derivatives 
are implemented using the central difference second order 
method - with second order accurate upwind/downwind 
at the boundaries. The first time derivative of A is needed 
to calculate J[b]i, and the first time derivative of J\b\iI^> 
appears in the Bondi news function. For both the quan- 
tities, apart from the first two time-steps, we use one 



sided- downwind second order accurate differentiation 
scheme. For the initial steps we use the standard first 
order scheme. 



V. COMPUTATIONAL TESTS AND RESULTS 

We tested our method for calculating the Bondi 
news function against two exact solutions: linearized 
Robinson- Trautman, and Schwarzschild in rotating co- 
ordinates. We discretize this analytic solutions on the 
required computational grid and determine the various 
metric components and use them as an input for the news 
module. The news module as such sees it as input from 
an equivalent numerical evolution code. 



A. Linearized Robinson- Trautman solution 

The Robinson- Trautman solution represents a dis- 
torted black hole emitting purely outgoing radiation. 
The radiation decays exponentially, and asymptotically 
the solution becomes Schwarzschild. 

The Robinson- Trautman metric can be put in the form 



ds 2 = -(JC-—)du 2 
rW 



2Wdudr - 2rW, A dudx A 

+r 2 q AB dx A dx B , (30) 



where K, = W 2 [l — L 2 (log>V)], L 2 is the angular mo- 
mentum operator (i.e., the angular part of the Laplacian 
operator on the unit sphere), and W(u,x A ) has to satis- 
fies the nonlinear differential equation 



12<9 n (logW) =W 2 L 2 IC 



(31) 

by putting 



We recover the Schwarzschild solution 
W ^constant in the above metric. 

Linearized solutions to the Robinson- Trautman equa- 
tion (|31|l are obtained by setting W = 1 + </> and dropping 
nonlinear terms in d>, 



I2d u 



L 2 {2-L 2 



(32) 



For a spherical harmonic perturbation <f> — 
X(u)Yi m this leads to the exponential decay 
A = \(0) e -^e(i+i)(i 2 +i-2)/i2^ The ( linear ) p ertur bation 

can be a linear sum of various spherical harmonics with 
small amplitudes. 

We consider a linearized solution of the form: 



= m*22e- 2u Y 22 + A 33 e- 



(33) 



with A 22 = 3 x 10 -7 and A 33 = 7 x 10~ 7 . The metric 
components are, J = 0, /3 = 0.50, U — d(j)/r (so U = at 
r = oo) and W = 1. Giving these quantities as inputs at 
the grid points at each time level, we call the news module 
to calculate the Bondi news function. We calculate the 
L 2 norm J[b]o (which should be zero), as well as the 
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Li norm of the error in the Bondi news function, whose 
analytic value is 



N(u,x A ) = - W _1 3 2 VV. 



(34) 



The computation was performed for the following grid 
sizes 



(a) 


Aq = 


Ap = 


(). 10(1 Ax = 


1/60, Au = 


-- 0.04 


(b) 


Aq = 


Ap = 


0.050, Ax = 


1/120, Au 


= 0.02 


(c) 


Aq = 


Ap = 


0.025, Ax = 


1/240, Au 


= 0.01 


(d) 


Aq = 


Ap 


0.020, Ax = 


1/300, Aii 


= 0.008(35) 



In this case as C[b] — C we ml d that ||«/[s]ol|2 constant 
with u and we plot | j«7[s]ol la at u = 1.6 for the discretiza- 
tions (l3~5|) : the observed convergence rate is 1.99. The 
Li norm of the error in the news is plotted against u 
for the discretizations (|35|l in Fig. We see convergence 
that, at earlier times, is approximately second order, but 
at later times saturation is observed. The convergence 
rate between different grids as u increases is given in Ta- 
ble |Q}. Further, Fig. [SJshows that the error for the finest 
grid [33 (d) saturates at about e~ 24 , the next finest grid 
OHJc) saturates at e -24 5 , grid l3*5T b') saturates at e -26 5 , 
while grid !35r a'l saturates at e -28 5 . We observe that the 
saturation level scales approximately with discretization 
length as A~ 3 . These matters are discussed further in 
Sec. ED 

The Y33 component initially dominates the total news 
function, which should therefore decay approximately as 
e -i0u j e w ith slope 10 on a log scale; while at later times 



convergence rate at u 


grids (a) & (b) 


(b) & (c) 


(c) & (d) 


0.4 


1.955 


1.992 


1.996 


0.8 


1.929 


1.982 


1.966 


1.2 


1.928 


1.966 


1.578 


1.6 


1.928 


1.907 


0.4746 



TABLE I: Change of convergence rate with u in the L2 norm 
of news function for linearized Robinson- Trautman for differ- 
ent grid resolutions in equation 1351 



the Y22 component dominates and the news should decay 
as e~ 2u . If the percentage of error remains constant in the 
news function, the Li norm of the error should also show 
a similar trend with increasing u. We see from Fig.^that 
this is actually the case. Until saturation effects become 
evident, a little before u — 2, the norm of the error in 
the news decreases as expected at all grid resolutions. 

The error in the news is significantly higher at the first 
two time-steps, presumably because we are using a first 
order scheme since there is not enough past data to eval- 
uate second order accurate time derivatives (actually, at 
the second time-step one can implement a second order 
scheme, but its error is larger). 



B. Schwarzschild in rotating coordinates 

By the transformation </? — > <p + nu of the azimuthal 
coordinate, the Schwarzschild line element in null coor- 
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-22 | 



_g g I i i i i i i i i 

' -4 -3.8 -3.6 -3.4 -3.2 -3 -2.8 -2.6 -2.4 -2.2 

ln(Aq) 

FIG. 3: Convergence of ||J[s]o||2 at it = 1.6 for the rotating 
Schwarzschild case. 

dinates can be written as 

ds 2 = -(1 - — - K 2 r 2 sin 2 6)du 2 - 2dudr 
r 

+2nr 2 sin 2 Odudip + r 2 q AB dx A dx B . (36) 

In this case we set n = 4 This coordinate change gives a 
nontrivial value for U at X + and thus it becomes a use- 
ful test to check that the numerically calculated Bondi 
news function remains zero. In this case the the values 
of the various metric components are /3 = 0, J = 0, 
U = 2IkQ/(\ + CC); W = 1. Giving these as inputs we 
calculate ||J[B]o||a and the news function ||iV||2, both of 
which have zero as the analytic value. We again found 
that ||J[b]o||2 hardly varies with u, and in Fig. |3]we plot 
ll^[B]ol|2 at u = 1.6 for the discretizations (j3*5)l : the ob- 
served convergence rate is again 1.99. The L2 norm of 
the news is plotted against u for the discretizations H35|) 
in Fig.Q] We find convergence to zero at u = 1.6 at a rate 
4.97 averaged over the grids (|35l) (with only minor vari- 
ations according to the specific grid used). The higher 
than expected rate is presumably due to some cancella- 
tion effect in the truncation errors. It is interesting that 
the magnitude of the error in the news builds up quite 
rapidly with time. Since the error remains convergent to 
zero, this must mean that the truncation error constants 
are growing with time. 

VI. CONCLUSIONS 

We have presented a method for calculating the news, 
based upon an explicit coordinate transformation from 



_2g 1 1 1 1 1 1 1 1 1 

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 

u 

FIG. 4: 1/2 norm of the error in the news as a function of time 
u for various discretizations for the rotating Schwarzschild 
case. 

general Bondi-Sachs coordinates to a Bondi-Sachs coor- 
dinate system that satisfies the Bondi conditions. The 
method has been implemented computationally and val- 
idated against known analytic solutions. There are sim- 
ilarities and differences between the method and that 
of 

• Both methods require J = on X + at u = 0. 

• Both methods solve evolution equations for Xr^ 
(Eq. d) and u [B]0 (Eq. JUS}). 

• The method of solves evolution equations for 
6 (Eq. (0-36)) and w (Given after Eq. (0-31)). 
In our method, S is not required, and u> is found 
explicitly (Eq. (JTSJl). 

• The formulas used here are much simpler than 
those of 0. 

• There is a means to monitor the accuracy 
of results obtained: we simply find J[b]o and 
compare the value to zero. 

The saturation and reduction in convergence rate re- 
ported in Sec. IV Al is typically an indication that round- 
off error becomes significant as the magnitude of the sig- 
nal decays. The saturation level scales approximately 
as A~ 3 , which corresponds to the finite difference rep- 
resentation of a third derivative. This is consistent with 
the fact that here the news calculation implicitly involves 
third derivatives: and U[b]o are differentiated to ob- 

tain A (Eg. I22f> ): then A is differentiated to obtain Jwu 
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(Eq. (125(1 ): and finally J\b\i is differentiated to obtain the 
news N (Eq. The observed scaling of saturation as 

A~ 3 is for the specific case of the Robinson- Trautman so- 
lution, which has the special property = x^; thus, 
in a general case the effect may be more severe. For future 
applications, it is important to be aware that round-off 
error may be a problem. 

It will be interesting to apply our news calculation to 
the problem considered in |l| (for which the method of Q 
was unsuccessful), but that is deferred to further work. 
Also deferred to further work is the possibility of extend- 
ing the coordinate transformation to higher order in £, 
and so being able to extract the Bondi mass and angu- 
lar momentum aspects - the difficulty will be that the 
computer algebra is likely to lead to very complicated 



formulas. 



Acknowledgements 

We benefited from the hospitality of the Max-Planck- 
Institut fur GravitationsPhysik, Albert-Einstein-Institut 
and the Center for Gravitational Wave Physics at Penn- 
sylvania State University (funded by the NSF coopera- 
tive agreement PHY-0114375). The work was supported 
by the National Research Foundation, South Africa un- 
der Grant number 2053724. We thank Sascha Husa and 
Jeffrey Winicour for discussions. 



[1] T.Regge and J. Wheeler, Phys. Rev. 108, 1063 (1957). 

[2] F. Zerilli, Phys. Rev. Lett. 24, 737 (1970). 

[3] V. Moncrief, Ann. Phys. (N.Y.) 88, 323 (1974). 

[4] M. Rupright, A. Abrahams, and L. Rezzolla, Phys. Rev. 

D 48, 044005 (1998). 
[5] K. Camarda and E. Seidel, Phys.Rev. D 59, 064019 

(1999). 

[6] E. Seidel and W.-M. Suen, Jnl. Comp. and Appl. Math. 

109, 493 (1999). 
[7] N. T. Bishop, R. Gomez, L. Lehner, and J. Winicour, 

Phys. Rev. D 54, 6153 (1996). 
[8] N. T. Bishop, R. Gomez, L. Lehner, M. Maharaj, and 

J. Winicour, Phys. Rev. D 56, 6298 (1997). 
[9] N. T. Bishop, R. Gomez, S. Husa, L. Lehner, , and 

J. Winicour, gr-qc/0301060 (2003). 



[10] Y. E. Zlochower, Ph.D. thesis, University of Pittsburgh 
(2002). 

[11] N. T. Bishop, Class. Quantum Grav. 10, 333 (1993). 
[12] N. T. Bishop, C. J. S. Clarke, and R. A. d'Inverno, Class. 

Quantum Grav. 7, L23 (1993). 
[13] R. A. Isaacson, J. S. Welling, and J. Winicour, J. Math. 

Phys. 24, 1824 (1983). 
[14] H. Bondi, M. J. G. van der Burg, and A. W. K. Metzner, 

Proc. R. Soc. London A269, 21 (1962). 
[15] R. K. Sachs, Proc. R. Soc. London A270, 103 (1962). 
[16] R. Gomez, L. Lehner, P. Papadopoulos, and J. Winicour, 

Class. Quantum Grav. 14, 977 (1997). 
[17] R. Gomez, Phys. Rev. D 64, 024007 (2001). 



